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ABSTRACT 

This report describes a new magnetohydrodynamic numerical model based 
on a hexagonal spherical geodesic grid. The model is designed to simulate astro- 
physical flows of partially ionized plasmas around a central compact object, such 
as a star or a planet with a magnetic field. The geodesic grid, produced by a 
recursive subdivision of a base platonic solid (an icosahedron) , is free from con- 
trol volume singularities inherent in spherical polar grids. Multiple populations 
of plasma and neutral particles, coupled via charge-exchange interactions, can 
be simulated simultaneously with this model. Our numerical scheme uses piece- 
wise linear reconstruction on a surface of a sphere in a local two-dimensional 
"Cartesian" frame. The code employs HLL-type approximate Riemann solvers 
and includes facilities to control the divergence of magnetic field and maintain 
pressure positivity. Several test solutions are discussed, including a problem of 
an interaction between the solar wind and the local interstellar medium, and a 
simulation of Earth's magnetosphere. 

Subject headings: magnetohydrodynamics (MHD) — methods: numerical — planets 
and satellites: magnetic fields — solar wind 
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Introduction 



Many astrophysical plasma processes occur is regions of space surrounding a central 
compact object such as a star or a planet. Examples include stellar winds, planetary 
magnetospheres, supernova blast waves, and mass accretion onto compact objects. In all 
these environments the central body (a star or a planet) is typically much smaller than the 
characteristic scales of the plasma flows. In solving this class of problem on a computer, 
radial grids are commonly used because resolution can be readily increased near the origin. 
The simplest and the most commonly used is the standard spherical polar (r, 9, ip ) grid 



[e.g., 



Washimi fc Tanaka 



1996 



Pogorelov fc Matsuda 



1998 



Ratkiewicz et al 



19981). This 



grid has a singularity on the z axis, where the control volume = r^Ar sin 6A6Aip, 
Ar, A^ and A(p being the grid cell dimensions in the radial, latitudinal, and azimuthal 
directions, respectively, vanishes as sin^ — > 0. For explicit methods, this requires a small 
global time step to satisfy the Courant stability condition for a system of h yperbolic 



conse rvation laws for the entire grid (implicit or semi- implicit methods (e.g., iToth et al. 



19981 ) don't suffer from this limitation). 



Spherical grids are also used in simulating the transport of energetic charged particle , 



20091). 



such as galactic cosmic rays, in turbulent astrophysical flows (iFlorinski fc Pogorelovl 
Transport models based on stochastic trajectory (Monte-Carlo) methods also suffer from 
the singularity on the polar axis. For example, in modeling cosmic-ray transport in the 
heliosphere it is common to align the z axis with the solar rotation axis. Because energetic 
particle transport (diffusion and drift) is very rapid at high latitudes due to a weaker 
magnetic field, a model must take vanishingly small time steps when a particle ventures 
close to the polar axis, which results in an inferior overall model efficiency. 

Time step requirements can be relaxed substantially by employing a grid that has a 
(nearly) uniform solid angle coverage. Examples include triangle, hexagon, or diamond 
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based geodesic grids (IDu et al. 



2003 



Yeh 



2007 



Upadhyaya et al 



2010f). obtained by a 



recursive division of a base platonic solid, and cubed sphere grids f lRonchi et al 



1996 



Putman Sz Lin 



20071 ). This paper introduces a framework for finite volume methods of 
solution of hyperbolic conservation laws in three dimensions, such as gas-dynamic or MHD 
systems, using spherical geodesic grids composed of hexagonal prism elements. Results are 
illustrated on a three-dimensional simulation of solar rotation and formation of corotating 
interaction regions (CIRs), an interaction between the solar wind and the surrounding 
local interstellar medium or LISM, and a simulation of Earth's magnetosphere. The new 
framework can be employed to model a broad range of large-scale 3D astrophysical plasma 



fiows around a compact object where high computational efficiency is a priority. 



2. Grid structure 

Our three-dimensional grid consist of a 2D geodesic unstructured grid on a sphere 
combined with a concentric nonuniform radial stepping with smaller cells near the 
origin. The 2D surface grid is a Voronoi tesselation of a sphere produced from a dual 
triangular (Delaunay) tesselation. The latter is generated by a recursive subdivision of 
an icosahedron. We use the geodesic grid generator software developed by the ICON 
project (h ttp://icon.enes.orq) ioT u se in atmospheric circulation modeling. An optimization 



algorithm ( iHeikes fc Randall 



19951), included in their code, produces a mesh with a 
difference in spherical surface areas between the largest and the smallest cells of less than 
10%. 

The number of vertices (V), edges (E) and faces (F) on a grid produced by /th division 
is given by 

iVv = 20 ■ 2^', iVE = 30 ■ 2^', iVp = 10 ■ 2^' + 2. (1) 



In this notation a level grid is dual to the original icosahedron projected onto a unit 




Fig. 1. — Prom left to right: level 3, 4, and 5 geodesic Voronoi grids produced by a recursive 
division of the base icosahedron. 
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sphere. The base shape of a control volume in a finite volume method is a spherical 
hexagonal prism with the exception of 12 pentagonal prisms located at the vertices of the 
original icosahedron. Level 3, 4, and 5 hexagonal geodesic grids are shown in Figured! 

A more detailed view of the Voronoi grid structure is given in Figure [2l A hexagonal 
face Fm (shaded) is shown surrounded by six adjacent faces Fmi-.F^e- The face centers 
of the dual triangle-based Delaunay grid (blue lines) are located at the Voronoi vertices 
(V); likewise, the former's vertices are at the Voronoi grid's face centers (F). The edges 
of the Voronoi and Delaunay grids on a sphere are mutually orthogonal and intersect at 
their midpoints. To achieve acceptable resolution in typical astrophysical fiow modeling 
problems level 5 or higher geodesic grids should be used. Most of our simulations use level 
6 mesh containing 40,962 Voronoi polygons. We emphasize that these 40,962 hexagons are 
distributed evenly over the surface of each spherical layer of cells. This gives us a resolution 
of about 3 X 10~^ steradians in solid angle which corresponds to an angular resolution of one 
degree. By condensing the mesh in the radial direction, one can achieve a further degree of 
refinement as required by the problem. 

We introduce a set of unit vectors normal to the edges of the Voronoi grid fimn, where 
the index m refers to the mth face and n = [1,6] is the number of the edge counted in a 
counter-clockwise direction. The unit vectors are tangential to the surface of the sphere. 
The corresponding edge lengths, measured along great circles, are designated Lmn- The 
outward radial unit vector at the cell center is r^- We also designate the area of the face 
on the unit sphere as A^- In this notation the control volume is equal to 

AVim = Amr^Ar^, (2) 

where i is the index on the radial axis. The areas Am are calculated by dividing each 
hexagonal face into six (five for pentagons) spherical triangles and adding up their areas 
using standard expressions from spherical trigonometry. Having defined our cell dimensions 



m4 



Fig. 2. — A close up view of the geodesic grid illustrating the relationships between a face 
Fm and its neighboring faces Fmi-.F^e- Unit vectors normal to the edges of the Voronoi face 
nmi..nm6 are shown. A fragment of the Delaunay grid is drawn with blue lines. 
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is this way we can proceed to integrate a system of conservation laws inside a control 
volume. 

It is worth pointin g out that a similar geodesic-mesh-based model was developed by 



Nakamizo et al. 



(120091 ) Their mesh was generated from a dodecahedron by first dividing 



each face into five triangles followed by a recursive subdivision of each triangle into four 



smaller triangles. The resulting unstructured grid top ology is similar 



Nakamizo et al 



but no t identical) to 
( 120091 ) computations 



our dual Delaunay grid. Interestingly, in the model of 
are also performed on a hexagonal grid, generated by connecting the centroids of the 
Delaunay triangles. 



3. MHD conservation laws 



For the heliospheric and magnetospheric problems (§6) we solve a modified set of MHD 
equations, written in terms of conservative variables U and fluxes F as 



where 



U 



pu 

e 



(3) 



v 



(4) 



pu 

puu + pl - BB/(47r) 
(e + p)u-B(u-B)/(47r 
uB - Bu 

in CGS units. Here p is density, u is velocity, B is magnetic field, I is a unit dyadic, 
P = Pg + is the total pressure, Pg being the gas kinetic pressure, and the energy 

density e is given by 

e = ^ + ^ + ^. (5) 



We employ two alternative models to control the divergence of magnetic field. The first 
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is the numerical scheme proposed by iPowell et al.l ( 119991 ) , where numerical magnetic field 
divergence is advected out of the simulation domain with the fiow velocity. This scheme 
modifies the system of conservation laws with a hyperbolic source term 

/ ^ 

B/(47r) 
u- B/(47r 
u 

The second scheme employs a generalized Lagrang e multiplier (GLM) for a mixed 



Q 



VB 



v 



(6) 



/ 



hyperbolic-parabolic correction (jPedner et al 
additional transport equation for ip 

dip 



20021 ) ■ The system (4) is extended with an 



dt 



+ V-(4B) = -^^, 



(7) 



where Ch is a constant, isotropic advection speed, taken to be somewhat fast er than the 



fastes t wave speed in the problem, and Cp is related to the rate of decay of ip. 



Dedner et al. 



( 120021 ) proposed two methods to fix the value of c^: (a) by fixing the time rate of decay 
of the GLM variable = exp(— Atc^/Cp), where At is the time step, and (b) by fixing 
the characteristic length over which the decay occurs, given by Id = c^/ch- Both methods 
are available in our code. In the GLM scheme the conservation law for magnetic field 
(Faraday's law) is modified to read 



(9B 

— + V ■ (uB - Bu + #) = 0. 

ot 



(8) 



The system (3) is integrated over a control volume AV^.^ shown in Figure |3] to obtain 
the finite volume method 



AKr 



At 



{r1^ii2^i+l/2,m — r1_ii2^ i-l/2,m) ' fm 
6 

'^i^'^i ^ ^ Lmn^ i(mn) ' Arrm ~l~ AVimQjm- 



(9) 



n=l 
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Fig. 3. — A prismatic control volume showing the unit vectors normal to the interfaces. 
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Here Fj(^„) is the flux at the center of the edge shared by the mth cell and its nth neighbor 
(where n= [1,6]). Note that the source term Q does not involve the right hand side of Eq. 
(7). The parabolic correction is applied by multiplying the value of ip obtained from the 
finite volume scheme (9) by the decay factor, 

rw, method (a), 

^|J^^|JX{ ' ^ ^' (10) 

^-Atch/id^ method (b). 

In our simulations we typically use 0.9 < < 1 and la equal to several times the smallest 
linear grid size. 

The divergence of the magnetic field is obtained from Gauss's theorem in the same way 
as V • F is calculated in Eq. (9). More generally, the divergence and curl operators acting 
on an arbitrary vector v may be written as 

/ 2 2 \ 6 -s 

_ ' (,'"j_|_i/2"^i+l/2,m ~ ^i-i/2^i-l/2,m) Lmn'^ran ' ^i{mn) ^ s 



n=l 



V X V = '- ^- '- h > ^ — -. (12) 

r-Ar- ^-^ r-A 

An evaluation of a curl is necessary when modeling energetic charged particle transport, 
where the particle's drift velocity is proportional to V x (B/B^). The values of primitive 
variables at face centers Vi±i/2 and Vi(m„) may be approximated as arithmetic averages of 
the values in the two cells separated by the face. A more accurate approach, adopted here, 
is to use interface resolved states obtained from a solution to the corresponding Riemann 
problem (see §5). 

The finite volume system of conservation laws is integrated with a second order unsplit 
TVD-hke method (see below). Right and left interface values are calculated in the usual 
way using some appropriate linear reconstruction to achieve second-order spatial accuracy. 
Fluxes are calculated from a solution to a one-dimensional (projected) Riemann problem 



at each cell interface. Finally, time is advanced using either a first order (Euler) or, more 
commonly, a second order (Runge-Kutta) scheme, depending on the nature of the problem. 



Reconstruction 



To achieve second order spatial accuracy we employ limited piecewise linear 
reconstruction on primitive variables V = (p, u, pg, B, t/^)"^. In the radial direction the 
simplest and the most robust limiter available is the MinMod, with slopes Sim given by 



= minmod(Si„, S 



(13) 



where the left and the right slopes on an asymmetric stencil are, respectively 

S,r = 2- 



l,m g-(- 2^*+li"^ '^i 



Ari_i + Ari 



Ari + Arj+i 



(14) 



Also available is the more compressive monotonized central (MC) limiter (Ivan Leer! 119771 ) 
with 



minmod 



(Ar, + Ar,+i)S-„ + (Ar,_i + Ar,)S+ 
Ari_i + 2Ari + An+i 

The third option is the weighted essentially non-oscillatory (WENO) limiter 



;WENO 



(15) 



(16) 



where the WENO weights w are given by 



(SL + e) 



where p is an integer constant here taken to be 2, and e is a small number, which we took 
to be 10~^^ in our simulations. The (S^)^^ and (S~„)^p terms are traditionally referred to 
as smoothness measures in WENO methodology. 

For two-dimensional reconstruction on the surface of a sphere the code can use either 



a minimum angle plane (MAPR) method (jChristov &: Popov 



20081 ) or a 2D version of the 
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weighted essentially non-oscillatory (WENO) scheme (iFriedrichI Il998l ) . Common to both 
schemes, a local two-dimensional coordinate system (^,?7) is introduced on the sphere, with 
its origin at the face center F^. The coordinates of the six adjacent cell centers {C,mn, Vmn) 
are then calculated in this frame. These coordinates are measured along two arbitrary great 
circles intersecting at right angles at the position of the central face Fm- The procedure is 
illustrated in Figure ID The angle A and the great circle distance between the face centers 
c are effectively polar coordinates on the surface of a sphere. The local coordinates of F^n 
are calculated as 

^mn = ccosA, rimn = csmA. (18) 

Next, the six two-dimensional slopes S^, S*^ are calculated from the triangles with 
vertices located at the cell centers F^, F^n, Fm,n+i (shown in blue in Figure [2]) as 

Vmn{^im,n+1 ^im) Vm,n+l{^imn ^im) ('\Q\ 

'lmn<,m,n+l <,mn'lm,n+l 

Sri ^mn{^ im,n+l im) ^m,n+l{^ imn ^im) fOC\\ 

imn — ~ 7 37 ~ • l^UJ 

'lmn^m,n+l 'imn'lm,n+l 

In the MAPR method we evaluated the average slopes as 



imn ' imn ' \ > 

The reconstructed slopes S^^^^''' are those given by Eqs. (19) and (20) for which the 
average slope (20) is the smallest. Thus the method is a two-dimensional equivalent of the 
MinMod limiter. In the WENO method the reconstructed slopes are weighted arithmetic 
averages of all six slopes, namely 

6 

qWENO^.r, _ - (nr)\ 

^im — 2^ ^imn^imn- K^'^j 
n=\ 

The weights Wj^n of each slope are calculated as 

^imn /ooN 

Wim„ = , (23) 

/ jr, ^imn 
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Fig. 4. — Calculation of adjacent cell center (F^„) coordinates in a local coordinate frame 
associated with a face F^. Every line shown is a segment of a great circle. 



- 15 - 



where 



+ e 



(24) 



where we again use p = 2 and e = 10~^^. Because edge centers he midway between the 
corresponding two face centers and F^n on the Voronoi grid, the reconstructed values 
at edge midpoints Vumn) can be computed as 



(25) 



These values are used to calculate the interceh fluxes according to Eqs. (4), (7), and (8). 

The code also implements a slope flattening algorithm that reduces the value of the 
slopes calculated by the reconstruction module in the vicinity of strong compressions 
(shocks). This prevents the occurrence of oscillations downstream of the shock. To 
construct a flattener, we calculate the minimum value of the fast magnetosonic wave speed 



in each computational cell and its neighbors in the same spherical layer and in the 



layers above and b elow (a total o : 



then calculated as ( Balsara et al. 



21 c ells). The shock detector function in each cell dim is 



20091 ) 



dir. 



mm 



1, 



(V ■ Vi)im^k 



mm X 
"'f,im" 



+ 1 



H 



(V ■ u)imAk 

mill X 



(26) 



where A/j^ is a characteristic dimension of the cell AVim, 5 is a constant of order 1 and 
H is the Heaviside step function. Subsequently, the slope in the cell im is calculated as a 
weighted sum of a slope obtained with the standard limiter, such as WENO, and that from 
a more diffusive limiter, such as MinMod or MAPR. For example, in the radial direction we 
could use 



WENO 



MM 



(27) 
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Riemann solvers 



The fluxes F are calculated from an (approximate) solution to the one-dimensional 
Riemann problem at each interface between the prismatic cells. A suitable solver must be 
able to handle supersonic and t ransonic flows without losing pos itivity. Our tests revealed 



that modern HLL-type solvers (IBatten et al. 



1997 



Gurski 



2004 ) were generally superior 



to other solver types for the solar wind-LISM interaction problem, where they were least 
likely to produce a negative pressure upstream of a very strong (Mach number > 10) 
shock. Genuinely m ulti-dimensional Riemann solvers are now appearing in the literature 



(IBalsara 



2010 



20121 ). and they offer substantial advantages on logically rectangular meshes. 
However, the analogous work for unstructured meshes is the topic of vigorous research and 
was not incorporated in the present work. 

An HLLC solver consists of four states: the left and the right unperturbed states plus 
two intermediate states separated by a tangential discontinuity. Designating the left and 
the right bounding wave speeds of the Riemann fan by Si and Sr, respectively, the intercell 
flux may be written as 

^ F,, Si > 0, 

F, + 5KUr-U0, Si<0<S\ 
F, + S',(U; -U,), S*<0<Sr, 

Fr, Sr < 0, 



(2^ 



where F; = F(Ui), F,. = F(Ur.), and S* is the speed of the intermediate wave (a tangential 
discontinuity). Because in a HLLC solver the normal velocity component and the total 
pressure only change across the outermost waves, the speed of the tangential discontinuity 
is readily calculated by applying the Rankine-Hugoniot conditions across these waves. This 
yields the speed 



S* 



PrUnA^r " Un,r) - Pr + Bl^J (47r) - plUn,l{Sl - Un,l) + Pi - Bl^J {Att) 
Pr{Sr - Un,r) - Pl{Sl - Un,l) 



(29) 
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where Un(i,r) and Bn{i^r) are the normal-projected velocities and magnetic fields in the left 
and right states, respectively. Suppose two prismatic cells and {1^1712) share an 

interface with an index rii = [1, 6] in the neighbor list of the first cell and ^2 = [1, 6] in the 
neighbor list of the second cell (the definitions of "first" and "second" are arbitrary here; 
they could be defined, for example, by using the condition that mi < 1712)- Then a normal 
velocity projection is defined as 

where Uj(mini) and Uj(m.2n2) are the reconstructed velocities given by (25), and n„^„^ is the 
unit vector normal to the interface rii of the cell {i,mi), pointing outward. The values for 
Bn(i,r) are computed in the same way. 

Several HLLC MHD solvers may be found in the literature, distinguished by their 
choice of the tangential velocity and magnetic field components in the intermediate states 
U;*^ (unlike in gas dynamic s, these states are not unique in MHD). Currently we employ a 



solver proposed by 



Lil ( 120051 ) ■ Its main feature is that no jump in magnetic field is permitted 



across the tangential discontinuity. 

I 

second option available in our model is the HLLD Riemann solver flMiyoshi fc Kusano 



20051 ). This type of solver incorporates two additional states U^** separated from the 
corresponding "single star" states by rotational (Alfvenic) discontinuities, propagating to 
the left and to the right of the middle wave with speeds SI and S* respectively, given by 

Si-S - (4^^*)i/2' \-S + (4^^*)i/2' l-^^i 
where 5* is given by Eq. (34) below, and 

5"; Un I ^ Sr 'Un,r /or>\ 

Pi = Pi g Pr = Pr g ■ (32) 
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The intercell flux is computed as 

' Fi, 



Fi + Si{V;-Ui), 
Fi + Si{VJ-Vi) + Sl{Vr-U 

Fr + Sr{v;-Vr) + s:{v*; 

Fr + SJV*-lJr), 



l)^ 



Si>Q, 
Si<Q< SI, 

SI <o< s*, 



s* <o< s;, 



S:<0< Sr, 



(33) 



S. < 0, 



The HLLD solver is somewhat less robust than the HLLC counterpart because it has a 
singularity when one of the extremal waves is a switch-on shock. When this condition is 
encountered, the program falls back to the HLLC algorithm which is singularity- free. 

Because of nonlinearity of the solvers, under rare circumstances one of the intermediate 
waves could fall outside of the bounding (fast) waves. To prevent this from happening, we 
take the speed of the bounding waves to be the maximum of the left, right, and intermediate 
HLL states. Because the HLL state depends on the wave speeds themselves, we perform 
an iteration procedure until the external waves could be moved out no further. In the 
event that either HLLC or HLLD solver fails to produce a positive pressure in any of the 
i ntermediate states, we fall back to the very robust but dissipative HLLE Riemann solver 



flEinfeldt et al. 



199ll ) with a single intermediate state U*. 



The methods described above are used without modification with the source term 
divergence cleaning algorithm (Eq. 6). However, the GLM method introduces two 
additional waves moving with the speeds ±c/j that carry changes in Bn and ip only. Because 
Ch is the fastest wave speed, these waves bound the "base" Riemann fan, comprised of 2, 3, 
or 5 waves in the HLLE, HLLC, and HLLD solvers, respectively. The intermediate states 
are readily obtained from the Rankine-Hugoniot conditions at the bounding waves as 

Bn,l + Bn,r A " 



B* 



2c, 



4.233 5.413 6.593 7.773 8.953 10,13 11.31 



Fig. 5. — Magnetic field magnitude at time t = 0.07 from the blast wave problem. 
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r 



(34) 



2 2 

These intermediate states serve as both the right and the left states for the actual Riemann 
solver. The advantage of this approach is that any possible jump in the normal component 
of B is taken up by these additional external waves. 

Very few genuinely three-dimensional test problems are available for spherical grids. 
For code verificatiqn we used a 3D blast w a ve pr oblem similar to those presented by 



Gardiner fc Stonel fl2008h and 



Balsara et al. 



(|2009[ ). The simulation region is constrained 



between rmin = 0.01 and r^ax = 0.5. The radial cell width Ar increased outward 
monotonically from 0.00052 to 0.047. The inner boundary was treated as a perfectly 
conducting sphere with reflecting boundary conditions imposed. The initial conditions are 
p = 1, u = 0, and p = 10 (r < 0.1), p = 0.1 (r > 0.1). The initial magnetic field is given 
by the standard potential solution for a perfectly conducting sphere in a uniform external 
field, namely 

^3 



Br = Bn 



B0 = -Bo{l + 



2r3 



cos^. 



sin^^. 



(35) 
(36) 



This solution was rotated such that the external field pointed in the direction 

(l/VS, 1/^3, 1/^3). We used Bo = 10 and 7 = 1.4. The system was evolved until t = 0.07. 

For this problem we chose a level 6 grid with 256 cells in the radial direction. The 
GLM version of the numerical scheme was used with the HLLC Riemann solver and 
WENO reconstruction. Figure |5] shows the magnitude of magnetic field at the end of the 



simulation on a linear scale. 



Gardiner fc Stonel (|2008|) and 



he flow struc t ure o f the solution is qualitatively similar to 



Balsara et al.l (120091 ). consisting of an outermost fast shock 



wave and two dense shells of material elongated along the magnetic field. This problem did 
not trigger the slope flattening or positivity correction routines meaning it is not a very 



good stress test of the code. Several more difficult problems simulating actual astropliysical 
plasma flows are discussed next. 



6. Numerical solutions of test problems 

To illustrate the capabilities of the new model we present results from three different 
simulations of solar system plasma environments. The first is a dynamic MHD simulation 
of compressive structures in the solar wind known as corotating interaction regions (CIRs). 
This is a simple test problem with a strong degree of spherical symmetry. The second is a 
simulation of the structure of the global heliosphere, including regions on each side of the 
interface between the solar wind and LISM known as the heliopause. This problem involves 
mode complex transonic fiows and a population of neutral atoms in addition to the plasma. 
Finally, our third test problem is a stationary structure of the Earth's magnetosphere. 
Unlike the two previous cases, this one is an example of a highly magnetized plasma 
environment. 



6.1. Test problem 1: Corotating interaction regions in the solar wind 



Corotating interaction regions (CIRs) are compressive structures produced through 
an interaction between high and low speed s treams in the solar wind. CIRs are fully 



formed by the time they reach Earth's orbit (jSiscoe 



1972 



Gosling et al. 



1972|). When the 



streams emanating from the Sun are approximately steady in the co-rotating frame, these 
compression regions form spirals in the solar equatorial plane that co-rotate with the Sun. 
The leading edge of a CIR is a forward compressional wave propagating into the slower 
solar wind ahead, whereas the tailing edge is a reverse wave propagating back into the 
trailing high speed stream. At large heliospheric distances the waves steepen into forward 
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and reverse shocks. The entire plasma structure is convected with the solar wind and plays 
an important role in the dynamics of the heliosphere. 



CIRs have been extensively studied usin g global MHD simulation ( jPizzol Il994 : 



Riley et al. 



2001 



Usmanov fc Goldstein 



20061 ) ■ To ge nerate CIRs in a global MHD 



simulation we adopt the tilted-dipole flow geometry of iPizzd (119821 ) at the inner boundary, 
which is illustrated in Figure |6l In this figure Oxyz is the fixed (heliographic) frame, where 
z is the solar rotation axis, and Ox'y'z' is a frame aligned with the Sun's magnetic axis z'. 
The parameter 7 is the dipole tilt angle, and /3 is the latitude of the fast-slow transition 
boundaries (blue circles) in the coordinate system Ox'y'z'. In the simulation discussed below 
we used /3 = ±30°. 



Following iPogorelov et al.l (120071 ) . one readily derives a quadratic equation for the 
latitude of the transition line ^ as a function of the azimuthal angle cp, 



asm 6^ + 6 sin 6^ + c = 0, 



(37) 



where 



a = cos^ 7 + sin^ 7 tan^ v9(cot 7 C0S7 + sin7)^ 

b = 2 sin /3 cos7(l + tan^ (y?) 

c = sin^ /3(1 + tan^ v9cos^7) — cos^ /3sin^7tan^ 



(38) 



Note that when /3 = 0, 9{ip) reduces to the express ion for the latitude of the magnetic 



equator given by Eq. (A6) of 



Pogorelov et al. 



(l2007f ). 



We simulated a region of the solar wind between r^in = 0.5 AU and rmax = 30 AU 
using 512 concentric grid layers of variable thickness (increasing outward). At 1 AU we 
assume the following conditions: density n = 3.5 cm~^ and radial velocity u = 800 km/s 



in the fast solar wind and n = 7 cm 



u 



400 km/s in the slow solar wind. The radial 



component of the magnetic field at 1 AU was Bf. = 28 fiG. These conditions were extended 
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Fig. 6.— 
tion. 



A diagram of the assumed titled-dipole plasma flow geometry for the CIR simula- 
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to the inner bo undary 



magnetic field (jParker 



using the conventional Parker solution for the solar wind and its 



19581 ). The dipole tilt angle was taken to be 7 = 20°. The boundary 
(shear layer) between the fast and the slow solar wind flows was located at a latitude of 
30° in the coordinate system aligned with the dipole axis. This simulation was performed 
on a level 6 geodesic grid. We chose the HLLC solver to evolve the time-dependent MHD 
equations, combined with the GLM divergence cleaning method; WENO reconstructions 
was used in all directions. 

Figure [7] (left) shows the logarithm of the magnetic field magnitude in the xz and xy 
planes using a cutout plot. Plasma velocity vectors are shown as arrows of variable length. 
The CIRs can be visually identified as higher density and magnetic field intensity regions 
(red). The maximum latitudinal extent of CIRs is given by the sum of the angle between 
the rotation and the dipole axes and the extent of the slow solar wind in the frame aligned 
with the dipole axis, i.e., 7 + /3 = 50°. In the equatorial plane, the spiral CIR structure is 
seen to be bounded by shock-like discontinuities. 

Several characteristic CIR features can be recognized in the plasma radial profiles 
shown in Figure [7] (right). We chose the profile along the direction 25° northern latitude 
relative to the solar equatorial (xy) plane. The forward-reverse sh ock pairs are common ly 



observed at mid-latitudes, below the heliographic latitude of 26° (iGosling &: Pizzo 



19991) 



They are shown with vertical dashed lines in the Figure. Shock pairs associated with 



CIRs are believed to 



3e responsible 



cosmic-ray intensity (iKota &: Jokipii 



' or the observed 2 6- day recur rent decreases in galactic 



1991 



McKibben et al 



19991 ). Other features, such 



as the south-north flows are also identified through the north-south flow deflection angle 
e = sin^^(— Me/lul) shown in the bottom panel. The transitions from northward (positive) 
to southward (negative) velocity are separated by roughly one Carrington rotation period 
(26 days) in our simulation. We conclude that the model is capable of reproducing the 
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5 10 15 20 25 

Heliocentric Distance (AU) 



Fig. 7. — Left: Magnetic field magnitude (log scale) in the meridional plane (xz) and the 
solar equatorial plane (xy) for the CIR simulation. Arrows are the plasma velocity vectors. 
Right: radial profiles along the direction {6, ip) = (25°, 135°) of (from top to bottom): 
plasma density (log scale), radial velocity, log thermal pressure, log magnetic field intensity, 
log temperature, and the north-south flow deflection angle e. Arrows mark the forward 
(pointing right) and reverse (pointing left) propagation of wave fronts. 
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essential CIR features and is consistent with the earher simulations of this phenomenon. 



6.2. Test problem 2: The global heliosphere 

The energy density in a supersonic stellar wind, such as the solar wind, decreases in 
inverse proportion to the square of the distance from the star. Eventually the outflow is 
unable to maintain pressure balance with the galactic environment near the star, comprised 
mostly of partially ionized hydrogen gas. The stellar wind undergoes a transition to a 
subsonic flow at a structure called a termination shock. A tangential discontinuity called 
an astropause (heliopause for the solar wind) separates the shocked stellar flow from the 
interstellar gas. A bow shock may develop in front of the astropause if the relative motion 
between the star and LISM is supersonic. In the case of heliosphere, the region between the 



termination shock and the boundary is called the heliosheath. The theor y of ste 
interfaces (as appli ed primarily to the he liosphere) has been developed in 



lar wind 



AxfordI (11972[ l. and 



Baranov et al. 



the interface could be found in 



Parked (119611 ) 



(119761). Recent thr ee-di mensional MHD sim ulations of 



Pogorelov et al. 



(I2007h and 



Opher et al. 



(I2007h . 



To simulate the structure of the heliospheric interface we used a relatively coarse level 
5 geodesic grid with 240 radial points. As in the CIR problem, the concentric layer spacing 
was nonuniform with the smallest cells at the inner radial boundary located at 10 AU; the 



outer boundary was placed at 900 AU. A heliogra phic coordinate syst em is used here, where 

me] 



the z axis is aligned with the solar rotation axis (iBeck &: Gilei 



20051 ). an d the x axis is in 



the p lane formed by the z axis and the interstellar helium flow direction ( iLallement et al. 



20051 ). The y axis completes the right-handed orthogonal system. The geometry of the 



problem is illustrated in Figure [HI 



The heliospheric configuration computed here is representative of a solar minimum 
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Fig. 8. — The heliographic coordinate system used in the simulation of the global heliosphere. 
The directions of the flow of interstellar hydrogen (Vh) and helium (Vhc) span the so-called 
hydrogen deflection plane (HDP) with a normal n. Here the interstellar magnetic field B 
lies in the HDP, with an angle of 45° relative to ^ He- 
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( lFlorinskill201ll ). At 1 AU we assume the following conditions: density n = 3.5 cm~^ and 
radial velocity u = 800 km/s at heliographic latitudes above 30° (fast solar wind) and n = 7 
cm~^, u = 400 km/s at low latitudes (slow solar wind). The magnetic field is a Parker spiral 
with a radial component Br = 28 at 1 AU. The azimuthal magnetic field component 
is a function of the solar wind speed. The heliospheric current sheet is not included in 
this simulation, so that the solar magnetic field is alw ays directed outward from the Sun. 



The observed current sheet is between 10^ km (1 AU 



times 10^ km (heliosheath, 



Burlaga fc Ness 



Winterhalter et al. 



19M) and a few 



201ll ) in width, which is much too narrow to be 



resolved with a global model. 

The interstellar flow has a total density of 0.2 cm~^, and is ionization rate of 0.25. 
Its velocity vector is Vhg = (—26.3, 0, —0.23) km/s in the chosen heliographic coordinate 
system. The interstellar magnetic field lies in the so-called hydrogen deflection plane (the 
plane spanned by the velocity vectors of neutral interstellar hydrogen and helium) and is 
inclined by 45° with respect to the LISM flow vector. Its components are (—1.3, 1.38, —2.32) 
yuG in our coordinate system. The temperature of both ionized and neutral components 
in the LISM is taken to b e 6530 K. Th e neutral and the plasma fluids are coupled via the 



charge exchange process (lAxford 



19721 ). We simulate both fluids using the same code by 



ex plicitly flxing B = for the neutral hydrogen. The charge exchange terms used are those 



of 



Pauls et al. 



( 119951 ). For simplicity we only include interstellar hydrogen in this simulation 
and ignore atoms produced by charge exchange in the heliosheath or the solar wind. To 
separate the interstellar region from the heliosphere we use a passively advected indicator 
variable q which satisfles the equation 

d{pq) 



dt 



+ V ■ (pgu) = 0. 



(39) 



The indicator variable is set to 1 in the solar wind and —1 in the interstellar flow. The 
condition g = then gives the location of the heliopause. 
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Fig. 9. — Left: Constant plasma pressure surfaces (log scale) cut by the meridional (xz) plane 
for the heliosphere simulation. Selected magnetic field lines in the LISM are shown. Right: 
Radial profiles in the upwind (nose) direction of (from top to bottom): log plasma number 
density, radial speed, log thermal pressure, magnetic field magnitude, and log temperature. 
The positions of the termination shock and the heliopause are marked with vertical dashed 
lines. 
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We chose the HLLC solver for this work because of its more robust handhng of a 
strong flow shear between the fast and the slow solar wind. We used the GLM V ■ B 
control method and WENO reconstruction in all directions. Simulations were run until a 
steady state was achieved which took about 300 years of simulated time. Figure [H left, 
shows a cutout view of the heliospheric interface. Surfaces of constant plasma pressure are 
plotted together with magnetic field lines in the LISM, illustrating their draping around the 
heliopause (the transition between the red and the green colors). The innermost pressure 
surface approximately traces the outline of the termination shock. 

We show radial profiles of several physical quantities in the upwind, or "nose" direction 
in the right panel of Figure |9l Before the termination shock, located at 67 AU in this 
simulation, the solar wind velocity is gradually decreasing because of a loss of momentum 
to charge exchange with interstellar hydrogen. In the heliosheath, the plasma density is 
nearly a constant while the magnetic pressure increases toward the heliopause where the 
flow becomes essentially stagnant. The effective heliosheath temperature (~ 3 x 10^ K) is 
that of the solar-wind and pic kup-ion mixture, which is significantly higher than that of the 



core solar wind ('~ 2 x 10^ K. 



Richardson et al. 



20081 ). From the top panel one can see that 



the density on the interstellar side of the heliopause is some 25 times higher than in the 
heliosheath. There is a very weak bow shock in this model barely visible in the pressure 
and temperature profiles. 

The results presented here were obtained using a single population of neutral hydrogen 
(the interstellar atoms). The computer code is actually capable of integrating conservation 
laws for multiple neutral hydrogen populations. It would be straightforward to include the 
heliosheath energetic neutral atoms and the neutral solar wind atoms in a simulation, at an 



added computational time expense (e.g. 



Williams et al. 



19971 ). 
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6.3. Test problem 3: Magnetosphere of Earth 



The Earth's magnetosphere is a product of an interaction between the supersonic 
solar wind and the geomagnetic field. Two major discontinuities, the bow shock and the 
magnetopause, are located between the undisturbed solar wind region and the geomagnetic 
field. The magnetosheath, filled with shocked solar wind plasma, lies between the bow 
shock and the magnetopause, which is the external boundary of the magnetosphere. The 
magnetopause thus separates the hot, tenuous magnetospheric plasma from the cold 
and dense solar wind plasma in the magnetosheath. Global MHD simulations, coupled 



with ionospheric mo dels, have been wide^ 
magnetosphere (e.g 



Fedder &: Lyon 



1995 



y used to study large-scale 



Tanaka 



1995 



Raedei 



1999 



processes in th e 



HuetaL 



20071 ). 



The geomagnetic field can be treated as a dipole field in the inner magnetosphere, its 
strength varying as r^^, where r is the distance from the center of the Earth. The thermal 
pressure varies more modestly leading to a very low plasma (3 (the ratio of the plasma 
thermal pressure to the magnetic field pressure) in the inner magnetosphere. Such low 
values of (~ ^Q"^ — 10~^) tend to produce numerical errors with conservative numerical 



schemes (IRaeder 



19991 ). To overcome this difficulty, the dipole field is treated a part from 



the to tal magnetic field according to the decomposition method introduced by 



Tanaka 



(119951 ). The momentum and energy fiuxes in the Riemann solvers are revised accordingly. 
The WENO reconstruction method is used in all directions and the GLM algorithm is used 
to control V ■ B. An interested reader will find more details on the GLM-MHD equations 
with a dipole field decomposition in the Appendix. 

The Geocentric Solar Magnetospheric (GSM) coordinate system is used in this 
simulation. It is centered at Earth, and the x, y, and z axes point to the Sun, the dawn-dusk 
direction, and along the north dipole axis, respectively. We choose the inner boundary 
to be a sphere with a radius r = 3Re (Earth radii), and apply the Dirichlet boundary 
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conditions. In particular, the number density is 370 cm~^, which is 1/27 of a typical value 
in the ionosphere. The thermal pressure is 4.65 x 10^^^ dyn/cm^, which is 9 times smaller 
than its ionospheric value. The magnetic field is taken to be a dipole field at the inner 



boun dary. For t 



(e.g., 



Janhunen 



le sa ke of simplicity, the magnetosphere-ionosphere electrostatic coupling 



19981 ) is not included, therefore the feedback of the ionosphere on the 
magnetosphere is ignored. We simply set the velocity to zero, which means there is no 
convection at the inner boundary. The free outer boundary is located at r = 100-R^;. 

We simulate a common configuration with a southward interplanetary magnetic field 
(IMF) of 50 fiG. The solar wind velocity is 600 km/s along the Sun-Earth line (the negative 
X direction), its number density is 5 cm~^ and temperature 9.1 x 10'^ K. The magnetic field 
is initially calculated as a superposition of a dipole field, centered at the origin, and a mirror 
dipole, located at (30-Re, 0, 0), The field on the sunward side is subsequently replaced with 
the solar wind field with = —50 fiG to make the initial configuration divergence free. In 
the simulation we used a level 6 geodesic grid and 256 grid points along the radial direction. 

A steady state configuration is obtained some 30 minutes (simulated time) into the 
simulation. The left panel of Figure [10] shows the color contours of the thermal pressure 
in the meridional (xz) plane and in the equatorial (xy) plane. The geomagnetic field and 
the IMF lines of force are also plotted. In the equatorial plane the geomagnetic field points 
northward, being opposite to the polarity of the southward IMF. We can see that the 
magnetosphere is open to the interplanetary medium and the geomagnetic field lines connect 



with the IMF ( iDungey 



19611 ). In that case the solar wind plasma momentum and energy 
can be transported into the magnetosphere through the site of magnetic reconnection. We 
did not observe surface waves or vortices induced by the Kelvin- Helmholtz instability (e.g.. 



GuoetaL 



20101 ) along the low-latitude magnetopause (the surface of the magnetopause is 



smooth in the equatorial plane) 
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Fig. 10. — Left: plasma pressure (color) and magnetic field lines for the magnetosphere 
simulation. Plane cuts for z = and y = are shown. The magnitude of the magnetic field 
is shown by the color of the field lines in the figure. Right: radial profiles along the Sun-Earth 
line of (from top to bottom): log plasma number density, velocity, log thermal pressure, log 
magnetic field intensity, and log temperature. The bow shock and the magnetopause are 
marked by vertical dashed lines. 
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The profiles of tlie pliysical quantities along the Sun-Earth line are shown in the 
right panel of Figure [TUJ The magnetopause is located at the neutral point for southward 
IMF case. The x velocity component approaches zero at the subsolar point, where the 
Sun-Earth line intersects the magnetopause. The shocked plasma becomes dense and hot 
in the magnetosheath, compared with the undisturbed solar wind. For southward IMF, the 
neutral point is found from the magnetic field strength profile (fourth panel from the top), 
where magnetic reconnection could occur in the presence of dissipation. 

Our result has all the relevant features of a typical MHD magnetospheric simulation. 
In this illustrative solution, we only calculate a steady state representative magnetosphere. 
Of course, the model can be also used with more realistic time-dependent IMF conditions 
derived from observations. 

7. Summary 

In this report we have presented a novel approach to numerical modeling of space 
plasma flows using geodesic spherical meshes with a nearly uniform solid angle coverage. 
This approach avoids the singularity on the symmetry axis inherent in polar spherical grids, 
leading to improved efficiency by allowing larger time steps. Our integration technique for 
gas-dynamic or MHD conservation laws is based on dimensionally unsplit time advance 
and uses two-dimensional reconstruction on the surface of a sphere. The new code has a 
number of useful features, such as a choice of multiple nonlinear Riemann solvers, weighted 
reconstruction limiters, and slope flattening to reduce possible oscillations near strong 
shocks. 

We have tested the new model on several common problems in space physics: a 
formation of corotating interaction regions in the solar wind, global modeling of the 
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heliospheric interface, and finally, the magnetosphere of a planet. Our results are consistent 
with those found in the literature and every feature of the resulting structures is well 
reproduced. At this time the model lacks an adaptive mesh refinement feature, which would 
permit a superior numerical resolution of shocks and discontinuities. Whereas a hexagonal 
(Voronoi) grid cannot be easily refined, its dual Delaunay grid can. The process starts 
with the original icosahedron that divides a sphere into 20 identical spherical triangles. 
Each triangle then may be recursively subdivided into four smaller triangles by connecting 
the midpoints of the original cell edges with great circle arcs. The Delaunay mesh is 
therefore naturally amenable to refinement based on an oct-tree formulation. Because each 
locally refined zone is further split in the radial direction, this is tantamount to each 3D 
patch giving rise to 8 identically-sized refined patches if it is to undergo one more level of 
refinement. 



The model could be potentially adapted to solve problems wher e the conipact 



Tanakal (|2000|), one 



object 



is not at the center of the region of interest. For example, following 
could introduce a non-concentric grid, where different spherical layer boundaries are offset 
from the origin. The offset distance increases for each subsequent layer, so that the mesh 
becomes denser in one direction and more rarefied in the opposite direction. Such an 
arrangement could be more efficient for modeling, e.g., a magnetosphere with a long tail. 

The new code by itself could be a valuable tool to investigate plasma flows around 
a source whose dimensions are small compared with the scale of the flow. Nevertheless, 
its chief intended purpose is to provide plasma background for subsequent simulations of 
the transport of energetic charged particles in the solar system and other astrophysical 
environments. Additional modules, recently added to the code, calculate the diffusion 
coefficients and drift velocity vectors based on magnetic field and other plasma properties. 
The use of geodesic grids will permit a more efficient calculation of phase space trajectories 
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in the stochastic integration method popular in cosmic-ray transport work (IBall et aLl 12005 



Florinski fc Pogorelov 



20091). The difference with polar grid-based models is expected to 
be quite pronounced in the polar regions of the heliosphere, where the diffusion and drift 
coefficients are typically very large. 

V.F. and X.G. were supported, in part, by NASA grants NNX10AE46G and 
NNX12AH44G, NSF grant AGS-0955700, and by a cooperative agreement with NASA 
Marshall Space Flight Center NNMllAAOlA. 

A. Dipole field decomposition 

In the magnetosphere, the external field Bi = B — B^, where B is the total magnetic 
field, and B^ is the internal dipole field. Since B^ is both curl-free (no current) and 
divergence free, we can write 

V • (b,B, - ^Bj?j = 0, (Al) 

(V X Bd) ■ (u X B) = 0. (A2) 
Using (Al) the momemtum flux from Eq. (4) can be expressed as 

puu + Pl - ^ (^BB + ^Bjl - B,B,^ . (A3) 

Next, from (A2) we obtain 

Brf ■ V X (u X B) - V ■ (u X B) X Bd = 0, (A4) 

which, upon substitution into the magnetic induction equation 

OB 

— -Vx{uxB) = (A5) 



yields 



We now define 



-V-[B(u-B,)-u(B-B,)] = 0. 



Bf pu^ pg Bl 

p* -p j^-L ei = - h-^ + — . 

Stt' 2 7-1 Svr 



Using these definitions the energy equation may be written as 

dei Brf dBi „ 



+— [u(Bi ■ Bi) + u(B . Bj) - B(u . Bl) - B(u • BJJ = 0. 

47r J 

Combining equations (A6) and (A8) yields 

^ + V ■ I (ei + pl)vi - ^[Bi(u ■ Bl) - u(Brf • Bi) + B,(u ■ Bi)] 



0. 



(A6) 



(A7) 



(A9) 



Using (A3) and (A9) the system of GLM-MHD equations with dipole field 
decomposition may be written as 



dt 



| + V.(p„) = 0. 
puu + p-I - -!- ( BB + \bII- BjBj 



(9Bi 



+ V ■ (uB - Bu + #) = 0, 



^ + V ■ I (ei + pt)u - ^ [(u ■ Bi)Bi - (Brf x u) x Bi] 



(AlO) 
(All) 
(A12) 
(A13) 

(A14) 



Note that the system (A10)-A(14) uses Bi and Ci instead of B and e as conserved quantities 



Consider the simplest three-state HLL solver (IHarten et al 
given by 

' F,, Si > 0, 

Fir, Si < < Sr, 
Fr, Sr < 0, 



19831). Its Riemann flux is 



(A15) 
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where Fi — F{\Ji) and = F(Ur) are the left and right unperturbed fluxes, respectively. 
The intermediate flux F;^ is given by 

SrFi — SiFr + SiSr(Ur ~ U/) 



■Ir 



(A16) 



Sr — Si 

Since only the definition of a conserved fiux is required to solve (A15), the system 
(A10)-(A14) can be readily used in place of (4). 

The decomposition of magnetic field does not affect the GLM scheme. For example, in 
the X direction we have two GLM equations. 

One can see that the external field component Bi^ can be integrated directly because the 
internal field (B,^ related terms) does not appear in these equations. 
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